Inclusion of the LithosphericTemp() option for the AddBox! function, e.g.

AddBox!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0),
        phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing), 
        DipAngle=30.0, T=LithosphericTemp(lbound="flux",nz=201))

The option enables to calculate a 1D temperature profile for a certain lithospheric setting (oceanic or continental) assuming variable thermal parameters (k,rho,cp) and radiogenic heating (H) for each phase. The temperature is then linearly interpolated from the 1D profile onto the entire box geometry, assuming a horizontal layering within the box. Rotation of the box is considered as well.

Therefore, one needs to define the phase distributions first in the AddBox! function using the Compute_Phase() function. The phases are defined as usual by phase = ConstantPhase() or phase = LithosphericPhases().

The thermal parameters for each phase are defined by the rheology structure, which is part of the LithosphericTemp() option.

The default values for the LithosphericTemp() option are:

    Tsurface  =     0.0             # surface temperature at the top of the box [C]
    Tpot         =     1350.0       # potential mantle temperature [C]
    dTadi       =     0.5             # adiabatic mantle gradient [K/km]
    ubound   =     "const"      # Upper thermal boundary condition ["const","flux"]
    lbound    =      "const"     # lower thermal boundary condition ["const","flux"]
    utbf         =      50.0e-3    # vertical heat flux q at the top [W/m^2]; if ubound = "flux"
    ltbf          =     10.0e-3    # vertial heat flux q at the bottom [W/m^2]; if lbound = "flux"
    age          =     120.0       # Lithospheric age for the temperature profile [Ma]
    dtfac       =     0.9           # Diffusion stability criterion for time iteration until age
    nz            =     201          # Vertical discretization of the 1D profile
    rheology =     example_CLrheology()      # Default thermal parameters assuming a three-layered lithosphere [Upper 
                                                                                         crust, lower crust, lithospheric mantle]

Following tests are included in test_setup_geometry.jl:
The figures show a model including a box, left: Phase distributions, middle: temperature field, right: Vertical temperature profile through the center of the box. The sum over T along this profile is also used for the test cases in test_setup_geometry.jl. If not otherwise defined, the test cases use the default values of LithosphericTemp().

1. a horizontally layer lithosphere; UpperCrust,LowerCrust,Mantle
   

2. an inclined lithosphere; UpperCrust,LowerCrust,Mantle
   

3. an inclined lithosphere with respect to the default origin (upper left corner of the box); UpperCrust,LowerCrust,Mantle
   

4. an inclined lithosphere with only two layers
   

5. using flux boundary conditions at the lower boundary of the box